Mathematical model combined with microdosimetric kinetic model for tumor volume calculation in stereotactic body radiation therapy

We proposed a new mathematical model that combines an ordinary differential equation (ODE) and microdosimetric kinetic model (MKM) to predict the tumor-cell lethal effect of Stereotactic body radiation therapy (SBRT) applied to non-small cell lung cancer (NSCLC). The tumor growth volume was calculated by the ODE in the multi-component mathematical model (MCM) for the cell lines NSCLC A549 and NCI-H460 (H460). The prescription doses 48 Gy/4 fr and 54 Gy/3 fr were used in the SBRT, and the effect of the SBRT on tumor cells was evaluated by the MKM. We also evaluated the effects of (1) linear quadratic model (LQM) and the MKM, (2) varying the ratio of active and quiescent tumors for the total tumor volume, and (3) the length of the dose-delivery time per fractionated dose (tinter) on the initial tumor volume. We used the ratio of the tumor volume at 1 day after the end of irradiation to the tumor volume before irradiation to define the radiation effectiveness value (REV). The combination of MKM and MCM significantly reduced REV at 48 Gy/4 fr compared to the combination of LQM and MCM. The ratio of active tumors and the prolonging of tinter affected the decrease in the REV for A549 and H460 cells. We evaluated the tumor volume considering a large fractionated dose and the dose-delivery time by combining the MKM with a mathematical model of tumor growth using an ODE in lung SBRT for NSCLC A549 and H460 cells.

various reactions in medicine, vaccine efficacy, and for predicting the effects of anticancer drugs [11][12][13] . An ordinary differential equation (ODE) in mathematical models is a set of differential equations containing an independent variable and one or more derivatives for that variable. An ODE is the most extensive form for modeling dynamical systems in science and engineering 10,14,15 . In systems biology, many biological processes (e.g., gene regulation and signal transduction) can be modeled by reaction rate equations that express the rate of the production of one species as a function of the concentration of another species in the system 11,13,14,16,17 . Evaluations of tumor growth using an ODE-based mathematical model in radiation therapy of tumors with single irradiation have been also been reported 18 . With the use of a mathematical model based on an ODE, it becomes possible to evaluate the tumor volume as an output, determine the effect of radiation therapy on the tumor, and derive the optimal irradiation schedule, irradiation dose, and the tumor effect when radiation therapy is combined with drugs [19][20][21] .
In order to compare the biological effects of radiotherapy administered with different doses and differing numbers of radiation treatments, mainly the following equation, called the linear quadratic model (LQM), is used in current radiotherapy 22,23 . Comparisons of the effects of the LQM have been made with single doses in the range commonly used in radiotherapy up to approx. 8 Gy 24,25 . In the LQM, the surviving fraction (SF) is calculated as a function of the absorbed dose (D) in Gy with two coefficients, α and β, where α is the proportionality factor to D [Gy −1 ] and β is the proportionality factor to D 2 [Gy −2 ]. However, the single fractionated dose exceeds 10 Gy in SBRT applied to lung cancers, and the actual cell survival rate may be higher than that predicted by the LQM when the dose is > 10 Gy [26][27][28][29] .
Some reports suggest that the LQM does not agree with the measured SF at high doses because the LQ curve bends continuously on a log-linear plot, which may interfere with extrapolation to high-dose fractionated treatments 28,30,31 . Various cell survival models have been proposed to solve this problem, microdosimetric kinetic model (MKM) was proposed that can predict the cell SF from physical doses based on domains, which are intracellular structures, for all types of radiation 32 . An MKM can accurately calculate cell viability even in the high-dose range by taking into account the radiation quality, the dose rate, and the cell DNA repair time [33][34][35][36][37] . There is a possibility of overestimating the effect of SF in the evaluation of LQMs, which may have a higher predictive at a single large dose since the calculation of cellular SF during irradiation by mathematical models based on ODEs is evaluated using LQM [19][20][21] . In addition, the LQM cannot account for sublethal damage repair (SLDR) affecting a tumor's survival during irradiation, and there are few reports of the effect of the dose-delivery time used to irradiate a prescribed dose on the tumor cell volume in mathematical models based on ODEs.
We therefore created a new mathematical model by combining an ODE and an MKM, and we describe the model below. Based on our validation of the model, we propose that this combination model can be used to predict the tumor cell lethal effect of SBRT administered to NSCLC.

Methods
The SF calculations for 6MV photon beams using the MKM. The cell nucleus is divided into hundreds of independent regions, which are called domains in an MKM 32 . Irradiation of these domains causes a potential lethal lesion (PLL). PLLs are classified into the following four categories according to their variants. (1) Irreparable lethal lesions (LLs) that appear in the primary process ('a' is the conversion rate constant); (2) PLLs that are converted to LLs in the secondary process ('b d ' is the conversion rate constant); (3) lesions that can be repaired in the primary process ('c' is the repair constant; and (4) lesions that do not become LLs for a certain period of time (t r ) and then become LLs and cannot be repaired.
The MKM assumes that a PLL is a double-strand break in DNA. The number of PLLs were caused a single instantaneous irradiation. The number of PLLs per domain using the rate constants of conversion (a, b d , c) for transformations is calculated as: Here, P is the number of PLLs in the domain and k d is the average number of PLLs per domain per dose [Gy −1 ] immediately after the irradiation. The parameter 'z' is the specific energy stored in the domain [Gy], and 't' is the time [h] after irradiation, satisfying 0 < t < t r . The t r is assumed to be infinite, as discussed by Hawkins 32 . The repair rate (a + c) of the PLLs was equal to the primary repair rate λ calculated by the DNA repair half-life. The average number of LLs (L n ) per cell nucleus ias defined as follows: While the corresponding parameters are.
Here, z is the specific energy [Gy] deposited in the domain, N is the number of domains, and A and B are coefficients. The parameter r d is the radius of the domain (0.5 μm), ρ is the density of the domain (1.0 g/cm 3 ), D is the absorbed dose (Gy), and y D is the dose mean lineal energy (keV/μm). The parameters α 0 and β 0 were determined by a single instantaneous irradiation using the LQM.
The MKM has been improved to take into account various dose rates and irradiation schemes with photon beams and changes in the amount of DNA per nucleus during irradiation 38 . The irradiation time and irradiation interruption time were considered for the cell SF, and the changes in the amount of DNA in the cell cycle were ignored (α 0 = constant, β 0 = constant); by taking the limit and setting N to infinity, the expression in Eq. (4) is transformed as follows: We thus transformed the Eq. (8) as follows: Here, Ḋ is the dose rate (Gy/min), and T is the dose-delivery time (min).
The parameter F is identified with the Lea-Catcheside time-factor G 39 . We used the α/β and DNA repair half-life T 1/2 values of the parameters for calculating the DNA repair rate (a + c) value of two NSCLC cell lines, A549 and NCI-H460 (H460) to evaluate the radiation effect of SBRT 28,29,40-42 . The dose mean lineal energy y D calculated by PHITS. The TrueBeam linear accelerator (Varian Medical Systems, Palo Alto, CA) using a 6MV X-ray beam was modeled with the particle and heavy ion transport code system (PHITS) 43 . PHITS can deal with photons, electrons, positrons, neutrons, and heavy ions 44 . The phase space files of the Monte Carlo (MC) that we applied were provided by Varian Medical Systems. The following phase space files were created using BEAMnrc built on the EGSnrc platform, and these phase space files created by BEAMnrc were transferred to the PHITS system in which the dose calculations were performed to simulate the 6MV photon beams. A water-equivalent phantom (20 × 20 × 20 cm 3 ) was created; the beam field size was 5 × 5 cm 2 with source to phantom surface distance (SSD) = 90 cm, and the measurement point was 10 cm deep. The calculation width was 3 cm in the water-equivalent phantom. The photon and electron cut-off energies were set to 0.01 MeV, and the MC calculation was performed with a statistical error < 1.0%. The dose-mean lineal energy y D for MKM was calculated as: where ε is the energy stored in the domain, l is the mean chord length, f(y) is the probability density of linear energy, and d(y) is the dose distribution of linear energy. The T-SED function of PHITS was used to calculate the dose mean lineal energy for 6 MV photon beams 45,46 . T-SED is a track structure T-SED calculates the distribution of energy imparted in a small area using formulas constructed based on the results of the analysis. The computationally derived y D values were used to evaluate the cell survival of the tumor in Eq. (4). www.nature.com/scientificreports/ Tumor growth volume calculation using the ODE with mathematical model. We used the ODE with a multi-component mathematical model (MCM) that models growth by distinguishing the state and growth rate of tumor cells (active tumors and quiescent tumors, etc.) in this study 20 . The MCM distinguishes the tumor cells into active tumors, resting cells that can stop dividing and turn into active cells, and non-dividing cells, which are dead cells waiting to be excreted into the bloodstream. The effects on tumor cells were evaluated by combining the growth of tumor cells represented by the MCM with the mortality rate of cells expressed using the radiation calculation. The tumor cells were divided into the types T 1 , T 2 , and …Tm as active tumors in the MCM. The volume of the T 1 tumor cells is V 1 ; the volume of quiescent cells (T Q ) is V Q , and the volume of the non-dividing cells (T ND ) is V ND . The radiation affects active tumors (T 1 , T 2 ,…T m ) and quiescent cells (T Q ), but not non-dividing cells (T ND ). Because cells in the same tumor might be in different states and have different growth rates, such as active tumors and quiescent cells, the following model can be constructed using Eq. (14).
The value of K 1 is constant related to the tumor growth rate and the environmental carrying capacity. The parameter p Q1 represents the probability that a T Q tumor cell is transformed into a T 1 cell, and p 1Q is the probability that a T 1 tumor cell is transformed into a T Q cell; p 1ND is the probability that a T 1 tumor cell is transformed into a non-dividing cell T ND in the blood. The growth of tumor cells represented by the MCM is modelled by Eq. (13). The volume of V 2 is evaluated in the same way using Eq. (14). The volumes of V m , V Q , and V ND are expressed by the following equations, since the active tumor is divided into M pieces using Eqs. (15)(16)(17).
We used the MCM optimized with M = 2 to simplify the evaluation using this model. The percentages of T 1 , T 2 , and T Q tumor cells in the total tumor volume were defined as 50%, 20%, and 30%, respectively, with good agreement for the measurement value. MCM biological parameters were optimized and matched to measurement values of A549 and NCI-H460 (H460) cells using MATLAB software with the SimBiology toolbox 47,48 .
The effect of SBRT on the tumor volume using the mathematical model combined with the MKM. The LQM has been used to calculate the cell lethal effect of photon beams [19][20][21] . In the present study, to more accurately assess the cell lethal effects of SBRT, we evaluated the cell lethal effects using the MKM instead of the LQM.
The lethal effects of radiation on the tumors were calculated in the model using the MKM, and Eq. (18) was converted to ODE format (Eq. (19)).
Here, Ḋ is the dose rate of radiation, and V is the volume of the tumor cells. The combination of tumor cell growth as indicated by the MCM and cell lethality as indicated by the radiological calculation yields the tumor cell volume calculation in SBRT (Eq. (20)).
Finally, to simulate fractionated irradiation in SBRT, the following equation was defined: www.nature.com/scientificreports/ V 1i , V Qi , and D i are the volume and absorbed dose of the i-th tumor cell at V 1 , V Q , and D, respectively. N is the total number of radiotherapy fractions. The prescribed doses in this study were 48 Gy/4 fr and 54 Gy/3 fr, based on the conditions of clinical trials of SBRT for lung cancer 49,50 . The tumor volume in the fractional irradiation of V 2 , V Q , and V ND was similarly calculated using Eq. (21). Figure 1 showed the MCM simulates combined with MKM to evaluate the effect of SBRT for tumor growth volume.
The dose rate Ḋ is 3 Gy/min, and V is the volume of the tumor cells, which represents the total tumor volume of V 1 , V 2 , and V Q combined, defined as 500, 1000, and 1500 cm 3 in this study. The t inter was set to 1 [day] in this study since t inter represents the time interval [days] of the fractionated dose (12 Gy or 18 Gy). The t intra represents the dose delivery time per one fractionated dose. Using the MKM instead of the LQM enabled the calculation of the effect of the dose rate to be taken into account.
In addition to the t intra derived from the relationship between the absorbed dose per fraction and the dose rate, we defined the t intra as 1, 5, 10, 30, and 60 min, and we then evaluated the effect of prolonging the dose delivery time t intra on the tumor volume. The first irradiation was defined as Time = 0 [day], the third irradiation, Figure 1. . We used the ratio of the tumor volume at 1 day after the end of the irradiation to the initial tumor volume to define the radiation effectiveness value (REV) in order to evaluate the effect of SBRT on the two NSCLC cell lines, A549 and H460.

Comparison of the measured tumor volume and calculated tumor volume by MCM in A549
and H460 NSCLC cells. Figure 2 shows the measured values of A549 and H460 and the tumor growth volume calculated using the MCM. The parameters of the MCM for tumor growth volume are given in Table 1.
The calculated results of the MCM for tumor volume evolution in time showed good agreement with the measurement values with both cell lines.
Comparison of the surviving fraction (SF) for the A549 and H460 cells predicted by the LQM and MKM models. The comparison of calculated SF by LQM and MKM for the measured tumor SF from irradiation of A549 and H460 cells were shown (Fig. 3). In the high-dose region (> 12 Gy), the LQM underestimate the SF compared to measured values, whereas the MKM calculations show good agreement in the highdose region. The parameters of the MKM for the tumor SF calculation are given in Table 2. The calculated cell viability with LQM and MKM for the measured cell viability with irradiation in A549 and H460 cells are shown (Fig. 3).   Table 3).

The impact of varying the V 1 , V 2 and V Q ratio on tumor volume in SBRT.
To determine its effect on the REVs in SBRT, we varied the ratios of V 1 , V 2 , and V Q tumor cells for the total tumor volume. Figures 5  and 6 depicts the results of the varying the ratio of V 1 and V 2 on the REVs. The percentage of V Q in the total tumor volume was fixed (30%, 500 cm 3 , 1000 cm 3 , and 2000 cm 3 ) and we varied the ratio of V 1 -V 2 at 10%, 20%, 35%, and 60% to evaluate the REVs.  3 . In addition, the tendency for an increase in the ratio of V 1 to increase the REVs were greater for the 54 Gy/3 fr compared to the 48 Gy/4 fr (Table 4). Figure 7a,b illustrates the effect on the REVs of changing the ratio of V Q . The REV in 48 Gy/4 fr of A549 with 10% V 1 (50 cm 3 ) 10% V 2 (50 cm 3 ) 80% V Q (400 cm 3 ) on 500 cm 3 initial tumor volume was 85.71% (Table 4).   www.nature.com/scientificreports/ The corresponding values with 20% V 1 (100 cm 3 ) 20% V 2 (100 cm 3 ) 60% V Q (300 cm 3 ) on 500 cm 3 initial tumor volume was 86.35% (Table 4). The higher the V Q ratio resulted in no significant effect on REV values. The effect of the V Q ratio on the REV was the same for the H460 cells and 54 Gy/3 fr. Figure 8a shows the effect of changing the dose delivery time per one fractionated dose t intra on the REV in A549 and H460 cells irradiated at 48 Gy/4 fr. The REVs of the A549 cells with 1 and 60 min t intra were 94.93% and 79.00% for 500 cm 3 , respectively (Table 5), and 95.75% and 85.46% for 2000 cm 3 , respectively. The REVs of the H460 cells with 1 and 60 min t intra were 96.39% and 88.80% for 500 cm 3 , respectively, and 96.91% and 91.59% for 2000 cm 3 , respectively. Figure 8b illustrates the effect of extending the dose delivery time per one fractionated dose t intra on the REVs for A549 and H460 cells irradiated at 54 Gy/3 fr. The REVs of the A549 cells with 1 and 60 min t intra were 99.15% and 92.37% for 500 cm 3 , respectively (Table 5), and 99.19% and 93.49% for 2000 cm 3 , respectively. The REVs of the H460 cells with 1 and 60 min t intra were 99.35% and 95.25% for 500 cm 3 , respectively, and   3 , respectively. The REVs were significantly decreased with the increase in the dose delivery time t intra ; in addition, the smaller the tumor volume, the greater the effect of increasing the t intra was. The effect of the extended irradiation time on the REV reduction was smaller for 54 Gy/3 fr compared to 48 Gy/4 fr. www.nature.com/scientificreports/

Discussion
In this study, we evaluated the effect of two different prescribed doses of SBRT on NSCLC cells by combining an MCM, which calculates tumor growth, and the LQM or the MKM, which calculate cell lethality from radiotherapy. As suggested in earlier studies [26][27][28][29][30][31] , there is a possibility of overestimating the SF in the high-dose range > 10 Gy when the LQM is used, suggesting the possibility of good SF prediction by using the MKM as reported (Fig. 3) 37,38,51,52 . In the present study therefore, the REVs for A549 and H460 cells had larger values for the LQM compared to the MKM: up to 4.9% and 2.7% higher at 48 Gy/4fr and 2.2% and 1.6% higher at 54 Gy/3fr www.nature.com/scientificreports/ (Table 3). The reason for the larger difference between the LQM and MKM REV in the A549 cells is thought to be that the difference in calculated SFs between the LQM and the MKM in the A549 cells was larger than that in the H460 cells (Fig. 3). In the case of 54 Gy/3 fr, the fractionated dose is larger than that of 48 Gy/4fr, which might explain the smaller difference in REV reduction between the LQM and MKM models, as shown in Table 3. A mathematical model combined with an MKM based on ordinary differential equations could thus be used to more accurately calculate the tumor volume for NSCLC in SBRT. Figures 5 and 6 showed the effect of varying the ratio of V 1 , V 2 on the REV. The difference in the REVs at 48 Gy/4 fr was up to 4.7% for A549 cells and up to 6.5% for H460 cells at an initial tumor volume of 500 cm 3 (Table 4), and the maximum was 3.1% for A549 cells and 4.0% for H460 cells at an initial tumor volume of 500 cm 3 at 54 Gy/3 fr. The greater the percentage of V 1 in the total tumor volume, the lower the REV was. The reason for this might be that V 1 represents the volume of active tumor T 1 , which has a smaller α/β value compared to active tumor T 2 (Table 1). In radiobiology, the effect is higher for a larger fractionated dose when the α/β ratio is small 53 .
The change in the percentage of quiescent cells T Q (V Q ) also affected the REV ( Table 4). The low radiosensitivity of quiescent cells compared to active tumor cells 54 . The advanced tumors have quiescent cells that grow slowly, and that the growth of human solid tumors depends not only on rapidly growing cancer cells, but also on their continued production 55 . Our comparison of the impact of V 1 and V 2 of active tumors and that of V Q of quiescent tumors as a ratio of the total tumor volume on the REV revealed that V 1 and V 2 had a greater impact on the REV. An estimated assessment of the ratio of active tumor to total tumor volume may therefore be needed for determining the REV; such an estimation is a task for a future study.
The flow cytometry was used to analyze the cell cycle distribution for the A549 and H460 lung cancer cell lines and observed that (1) the active tumors are highly dependent on the cell cycle, and (2) the rate of active tumors varies between tumor cell lines 56 . In addition, the tumor α/β values for the surviving fraction have been reported to vary even for the same type of lung cancer 57 . In the future, in order to optimize the parameters of the model, it will be necessary to collect experimental data and clinical data and optimize parameters such as V 1 , V 2 , V Q , and α/β.
There are several reports regarding the effect of a prolonged dose-delivery time (t intra ) on tumor cell survival, and a prolonged dose-delivery time was reported to decrease the effect on tumors 37,40,41,58 . The delivered biologically effective dose (BED) levels were calculated due to alterations in the SLDR, and a loss of BED was observed on 13 Gy with the dose-delivery time of 35 min compared to the acute-exposure approach. The effect www.nature.com/scientificreports/ of a prolonged dose delivery time using the MKM on tumors showed an approximate 6% decrease in relative biological effectiveness at an irradiation time of 60 min 40 . In the present study, a tumor volume of 500 cc at 48 Gy/4 fr and an irradiation time of 60 min showed the greatest REV reduction at 15.9% in A549 cells (Table 5).
On the other hand, H460 cells showed a 7.6% REV reduction under the same conditions. The greater impact of a time extension on the REV for A549 cells compared to H460 cells might be due to the faster DNA repair time for A549 cells and the greater impact of SLDR (Table 2). We suspect that the reason for the smaller REV reductions in both A549 and H460 cells at 54 Gy/3 fr compared to 48 Gy/4 fr is that the cell lethal effect of the larger single fractionated dose is less affected by the SLDR effect of the longer dose-delivery time. Considering the results of the single irradiation in other studies that evaluated a prolonged dose-delivery time together with the results of the fractionated dose used in the present study, our results showed a decrease in the REV, similar to the previous studies. Few investigations have taken into account the dose-delivery time of each fraction in a fractionated dose and calculated the effect of the fractionated dose on the tumor. Our present findings suggest that (1) the dose-delivery time might have an effect on tumors even in a fractionated dose, and (2) it might be necessary to attempt to increase the dose-delivery time as quickly as possible.
This study has two limitations to address. The first limitation is our evaluation of the REV by SBRT using ODE to estimate tumor growth and the cell lethality calculation model MKM, without considering changes in the tumor cell environment after irradiation. The high single radiation doses delivered, such as SBRT, may induce elevated and possibly persistent tumor hypoxia in NSCLC tumors 59 . Changes in tumor oxygenation after irradiation may alter the tumor response to radiation therapy (may impact effectiveness of radiation therapy); this requires further investigation 60,61 . The second study limitation is that we set 1 day after irradiation as the time point for the evaluation of the REV. For accurate REV derivation, it will be necessary to evaluate the effect on tumors with using more time evaluation points, as in clinical trials 8,46,47 .

Conclusions
We evaluated the tumor volume considering a large fractionated dose and the dose-delivery time by combining the MKM with a mathematical model of tumor growth using ODEs in lung SBRT for non-small cell lung cancer. Our results demonstrated that the ratio of active tumor to the total tumor volume and the dose-delivery time both affect the tumor volume in SBRT. www.nature.com/scientificreports/

Data availability
The data that support the findings of this study are available from corresponding author but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of corresponding author.  www.nature.com/scientificreports/